查看原文
其他

HiC数据分析实战之Hic-pro

生信技能树 生信技能树 2022-06-07

通过第一讲:三维基因组学习笔记,我们了解了3D基因组研究范围,然后根据我在生信技能树发布的生信工程师标准提炼出基础技能,也就是第二讲:生信基础技能 。最后提炼出了数据分析流程,并且安装好了对应的软件,也就是第3讲:流程及软件 。不过中间我还插播了一个文献解读 。

现在才正式开始数据处理实战,其中实战的测试数据,参考基因组以及对应的软件安装都是在第3讲:流程及软件 。看懂了这些准备工作,现在就可以跟我一起来一步步走通Hic数据分析流程啦,首先回忆一下准备工作咯。

conda create -n hic  python=2 bowtie2
source activate hic
conda search hiclab
conda install -y sra-tools samtools 
conda install -y bedtools trim-galore
source activate hic
conda install -y numpy scipy matplotlib h5py cython numexpr statsmodels  scikit-learn pandas 
cd ~ 
ln -s /data/training/zengjianming/ data
mkdir -p /home/zengjianming/data/biosoft/hiclib
cd /home/zengjianming/data/biosoft/hiclib
pip install https://bitbucket.org/mirnylab/mirnylib/get/tip.tar.gz
pip install https://bitbucket.org/mirnylab/hiclib/get/tip.tar.gz ##  17.7MB  , 44kB/s
mkdir -p /home/zengjianming/data/biosoft/hicpro
cd /home/zengjianming/data/biosoft/hicpro 
git clone https://github.com/nservant/HiC-Pro.git ##  46.95 MiB | 134.00 KiB/s, done
cd HiC-Pro/
mkdir /home/zengjianming/data/biosoft/hicpro/bin
mkdir -p /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0
mkdir -p /home/zengjianming/bin

再次了解 Hic-pro 软件

Hic-pro安装比较麻烦,我就不重复写了,大家仔细查看第3讲:流程及软件 。

很明显,前面教程我们讲到我们的Hic-pro安装在:

source activate hic
/home/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/bin/HiC-Pro -h

但是使用的学员服务器没有配置好,导致我的home目录下面空间不足,所以我挪动了data盘下面的空间,这样就有28T的空间啦。

所以软件挪了地方:

/data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/bin/HiC-Pro  --help

然后我们这一讲的重点是Hic-pro软件的使用。

如果仅仅是看软件的帮助文档,会发现其实非常容易使用;

usage : HiC-Pro -i INPUT -o OUTPUT -c CONFIG [-s ANALYSIS_STEP] [-p] [-h] [-v]
   -i|--input INPUT : input data folder; Must contains a folder per sample with input files
   -o|--output OUTPUT : output folder
   -c|--conf CONFIG : configuration file for Hi-C processing

重点就是3个参数而已,而最重要的当然是配置文件的准备咯。作者给出的示例文件不是一般的复杂;

## 记住  /data/training/zengjianming/ 目录就是 /home/zengjianming/data 目录。
cat /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/config-hicpro.txt
mkdir -p ~/data/project/hic/
cd ~/data/project/hic/

参考基因组

mkdir -p ~/data/project/hic/ref
cd ~/data/project/hic/ref
wget ftp://ftp.ensemblgenomes.org/pub/bacteria/release-40/fasta/bacteria_20_collection/caulobacter_crescentus_na1000/dna/Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa.gz
gunzip Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa.gz
source activate hic
bowtie2-build  Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa   bacteria
cat >genome.size
Chromosome 4016942

测试数据

来自于Tung B. K. Le et al. Science 2013 :https://www.ncbi.nlm.nih.gov/sra/?term=srr824846

mkdir -p ~/data/project/hic/fq/s1/
cd ~/data/project/hic/fq/s1/
858M Jul  3 16:21 SRR824846_Q20L10_1.fastq.gz
857M Jul  3 16:22 SRR824846_Q20L10_2.fastq.gz

运行Hic-pro软件

上面说过使用这个软件需要修改一个配置文件, 简单的包括:

BOWTIE2_IDX_PATH = # bowtie2建立的索引所在的路径,记住绝对路径
REFERENCE_GENOME = # bowtie2建立的索引
GENOME_SIZE = # 一个文件记录着参考基因组中每条序列的大小,格式为chr01 10000000

复杂的是GENOME_FRAGMENT =,也就是 HiC消化片段位点文件,这个文件在每个HIC流程中都需要生成,一般HiC建库的酶为hindiii(A^AGCTT ) 或者dnpiii 生成命令为 /PATH/HiC-Pro-master/bin/utils/digest_genome.py -r hindiii -o Refgenome.fasta

在软件例子里面是:

GENOME_FRAGMENT = HindIII_resfrag_hg19.bed
LIGATION_SITE = AAGCTAGCTT

但是我们这里的测试数据是 新月柄杆菌 Caulobacter crescentus , 文章里面 Le TB et al., "High-resolution mapping of the spatial organization of a bacterial chromosome.", Science, 2013 Oct 24;342(6159):731-4 提到的是限制性酶切位点NcoI 也是 用的内切酶,其它还有 BamHI、EcoRI、HindIII、NdeI、XhoI, 参考:https://www.biomart.cn/experiment/430/457/741/43956.htm

那么我需要针对特定酶和特定的参考基因组序列来设置好GENOME_FRAGMENT参数所需文件。

代码如下:

mkdir -p ~/data/project/hic/digest
cd ~/data/project/hic/digest
bin=/data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/bin/utils/digest_genome.py  
$bin -r C^CATGG  -o  bacteria.bed ../ref/Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa 

得到的 bacteria.bed 文件如下;

Chromosome    0   3188    HIC_Chromosome_1    0   +
Chromosome    3188    3551    HIC_Chromosome_2    0   +
Chromosome    3551    3737    HIC_Chromosome_3    0   +
Chromosome    3737    4500    HIC_Chromosome_4    0   +
Chromosome    4500    5578    HIC_Chromosome_5    0   +
Chromosome    5578    10834   HIC_Chromosome_6    0   +
Chromosome    10834   10840   HIC_Chromosome_7    0   +
Chromosome    10840   19183   HIC_Chromosome_8    0   +
Chromosome    19183   19381   HIC_Chromosome_9    0   +
Chromosome    19381   21967   HIC_Chromosome_10   0   +

然后修改配置文件后即可成功运行hicpro软件啦。

cd ~/data/project/hic/ 
source activate hic
cp /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/config-hicpro.txt ./
bin=/data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/bin/HiC-Pro
nohup $bin -i fq -o out -c config-hicpro.txt 1> run.log 2>&1 &

这个配置文件问题简单,重要的地方如下:

BOWTIE2_IDX_PATH = /home/zengjianming/data/project/hic/ref/ 
REFERENCE_GENOME = bacteria
GENOME_SIZE = /home/zengjianming/data/project/hic/ref/genome.size
PAIR1_EXT = SRR824846_Q20L10_1
PAIR2_EXT = SRR824846_Q20L10_2
GENOME_FRAGMENT = /home/zengjianming/data/project/hic/digest/bacteria.bed
LIGATION_SITE = CCATGG

但是 我遇到的报错居然是在 input的fastq文件的folder这个参数,太尴尬:https://groups.google.com/forum/#!msg/hic-pro/PrX00iA0N1U/apzlw8ArBwAJ 也就是说每个样本的fq文件都需要在自己独立的文件夹下面。

因为软件包含的步骤较多,所以会比较耗时!

输出结果解读

所有的输出,都在自己运行软件的时候指定的out文件夹下面。

共有bowtie_resultshic_results两个结果文件夹。

其中bowtie_results文件夹下共有三个文件夹:bwt2、bwt2_global和bwt2_local,分别是序列比对结果、染色体间关联比对结果和染色体内部关联比对结果。其中bwt2文件夹下有一些数据统计结果的输出文件,如mpairstat文件(如下)。

Total_pairs_processed    23546961    100.0
Unmapped_pairs    22197   0.094
Low_qual_pairs    0   0.0
Unique_paired_alignments    17744193    75.357
Multiple_pairs_alignments    462162  1.963
Pairs_with_singleton    5318409 22.586
Low_qual_singleton    0   0.0
Unique_singleton_alignments    0   0.0
Multiple_singleton_alignments    0   0.0
Reported_pairs    17744193    75.357

hic_results文件夹下共有data、matrix以及pic三个文件夹。

  • data文件夹下是比对上的有效序列对,文件末尾为_allValidPairs;

  • pic文件夹下是各类结果数据统计;

  • matrix文件夹下分为 iced 和 raw 两个文件夹,分别标准化后的关联矩阵和初始的关联矩阵。

重点必须是hic_results文件,如下:

hic_results/
|-- [  24]  data
|   `-- [4.0K]  s1
|       |-- [504M]  _bacteria.bwt2pairs.DEPairs
|       |-- [ 93K]  _bacteria.bwt2pairs.DumpPairs
|       |-- [ 83M]  _bacteria.bwt2pairs.REPairs
|       |-- [ 303]  _bacteria.bwt2pairs.RSstat
|       |-- [ 73M]  _bacteria.bwt2pairs.SCPairs
|       |-- [   0]  _bacteria.bwt2pairs.SinglePairs
|       |-- [1.2G]  _bacteria.bwt2pairs.validPairs
|       |-- [1.2G]  s1_allValidPairs
|       |-- [ 150]  s1_allValidPairs.mergestat
|       `-- [ 486]  s1.mRSstat
|-- [  24]  matrix
|   `-- [  41]  s1
|       |-- [  99]  iced
|       |   |-- [  85]  1000000
|       |   |   |-- [ 275]  s1_1000000_iced.matrix
|       |   |   `-- [ 125]  s1_1000000_iced.matrix.biases
|       |   |-- [  83]  150000
|       |   |   |-- [6.6K]  s1_150000_iced.matrix
|       |   |   `-- [ 675]  s1_150000_iced.matrix.biases
|       |   |-- [  81]  20000
|       |   |   |-- [345K]  s1_20000_iced.matrix
|       |   |   `-- [4.9K]  s1_20000_iced.matrix.biases
|       |   |-- [  81]  40000
|       |   |   |-- [ 86K]  s1_40000_iced.matrix
|       |   |   `-- [2.4K]  s1_40000_iced.matrix.biases
|       |   `-- [  83]  500000
|       |       |-- [ 792]  s1_500000_iced.matrix
|       |       `-- [ 225]  s1_500000_iced.matrix.biases
|       `-- [  99]  raw
|           |-- [  99]  1000000
|           |   |-- [ 139]  s1_1000000_abs.bed
|           |   |-- [ 165]  s1_1000000.matrix
|           |   `-- [  18]  s1_1000000_ord.bed -> s1_1000000_abs.bed
|           |-- [  96]  150000
|           |   |-- [ 783]  s1_150000_abs.bed
|           |   |-- [4.1K]  s1_150000.matrix
|           |   `-- [  17]  s1_150000_ord.bed -> s1_150000_abs.bed
|           |-- [  93]  20000
|           |   |-- [5.9K]  s1_20000_abs.bed
|           |   |-- [216K]  s1_20000.matrix
|           |   `-- [  16]  s1_20000_ord.bed -> s1_20000_abs.bed
|           |-- [  93]  40000
|           |   |-- [2.9K]  s1_40000_abs.bed
|           |   |-- [ 53K]  s1_40000.matrix
|           |   `-- [  16]  s1_40000_ord.bed -> s1_40000_abs.bed
|           `-- [  96]  500000
|               |-- [ 253]  s1_500000_abs.bed
|               |-- [ 469]  s1_500000.matrix
|               `-- [  17]  s1_500000_ord.bed -> s1_500000_abs.bed
`-- [  24]  pic
    `-- [ 188]  s1
        |-- [4.9K]  plotHiCContactRanges_s1.pdf
        |-- [5.2K]  plotHiCFragment_s1.pdf
        |-- [6.3K]  plotHiCFragmentSize_s1.pdf
        |-- [5.0K]  plotMappingPairing_s1.pdf
        `-- [5.1K]  plotMapping_s1.pdf

18 directories, 40 files

因为这个菌参考基因组很小,才4M左右,所以取500K,就只有8个bin了,其实意义不大。

随便秀一个pdf图吧。


似乎是还停留在对比对结果的统计层面上!


后续再更新结果的生物学意义解读,以及后续的compartment以及TAD和loop的分析,敬请期待哦!

里面提到的各种文献,综述,测试数据,软件,待这个三维基因组更新完毕会统一打包发放哦,请大家不要着急。

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存